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Venus Gravity Study Final Report 


This is a final report and the July 1991 semiannual report of the work 
done under NASA Grant NAGW-2361 for the preparation and publication of a 
report on the previous work that yielded a model of the gravity field of Venus 
based on the Doppler tracking data from the Pioneer Venus Orbiter. The grant 
for $5,960 was in response to our letter proposal P2340- 11-90 of 21 November 
1990. 


A paper describing our model of the gravity of Venus has appeared in 
JGR (Planets), copy appended. In order to make our work more useful to the 
community, we provided our digital map to the NSSDC, as indicated in the last 
paragraph of the paper. 

The paper in JGR appears to have been well received. In March 1993, I 
received a request from Dr. Peter Cattermole, FRAS, for permission to include 
the color maps of Venus gravity and smoothed topography in a forthcoming 
book, VENUS — THE NEW GEOLOGY, University College London Press. I have 
agreed to his request and sent him all necessary materials. 
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High-Resolution Gravity Model of Venus 

R. D. Reasenberg and Z. M. Goldberg 


Radio and Geoastronomy Division, Smithsonian Astrophysical Observatory, Harvard-Smithsonian Center for Astrophysics, 

Cambridge, Massachusetts 


The anomalous gravity field of Venus shows high correlation with surface features revealed by 
radar. We extract gravity models from the Doppler tracking data from the Pioneer Venus Orbiter 
(PVO) by means of a two-step process. In the first step, we solve the nonlinear spacecraft state 
estimation problem using a Kalman filter-smoother. The Kalman filter has been evaluated through 
simulations. This evaluation and some unusual features of the filter are discussed. In the second step, 
we perform a geophysical inversion using a linear Bayesian estimator. To allow an unbiased 
comparison between gravity and topography, we use a simulation technique to smooth and distort the 
radar topographic data so as to yield maps having the same characteristics as our gravity maps. The 
maps presented cover 2/3 of the surface of Venus and display the strong topography-gravity 
correlation previously reported. The topography-gravity scatter plots show two distinct trends. 


1. Introduction 

The Doppler tracking data from the Pioneer Venus Orbiter 
(PVO) contain the signature of the only available direct 
probe of the planet’s internal structure: gravity. Illuminating 
that internal structure is the prime objective of our study of 
Venus gravity. Early investigations of the gravity field of 
Venus were based on the analyses of tracking data from 
fly-by missions [Anderson and Efron, 1969; Howard et al., 
1974; Akim et al., 1978]. With the possible exception of the 
planetary mass estimate, these studies are superseded by the 
analyses of PVO data. In addition to the studies of local 
features by Phillips et al. [1979], Sjogren et al. [1980, 1983, 
1984], Reasenberg et al. [1981, 1982], and Goldberg and 
Reasenberg [1985], there have been evaluations of the 
low-order spherical harmonic components of the field by 
Ananda et al. [1980], Williams et al. [1983], and Mottinger et 
al. [1985] and a determination of an eighteenth degree and 
order spherical harmonic model by Bills et al. [1987], The 
gravity model presented here provides the highest resolution 
yet available, and the companion topography model, which 
mimics the resolution of the gravity model, provides a 
convenient means of forming and checking geophysical 
hypotheses. 

For an overview of the PVO mission, see Colin [1980]; for 
a further description of the scientific results, see the refer- 
ences and bibliography therein and the other papers in the 
same issue of the Journal of Geophysical Research. More 
recent results are given by Hunten et al. [1983]. 

Our analysis is based on a two-stage procedure. In the first 
stage, we use a Kalman-Bucy filter-smoother to solve the 
nonlinear spacecraft state estimation problem. The filter is 
incorporated in the Planetary Ephemeris Program (PEP) 
[Ash, 1972]. Significant attributes of this filter are described 
in section 4. The results of some studies of the filter’s 
performance are presented in section 5. Our analysis of the 
Venus gravity data is the first successful use of a Kalman 
filter as a tool for determining a planetary gravity field. 

The geophysical inversion performed in the second stage 
is described by Goldberg and Reasenberg [1985], who also 
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briefly discuss the Kalman-Bucy filtering. The inversion is 
carried out in a series of overlapping patches, each covering 
about 50° of longitude and extending as far poleward as the 
data allow. The overlap is sufficient to yield a band of at least 
10° longitudinal width in which there is no significant differ- 
ence in the models. Outside of that band, edge effects are 
noticeable, and the solutions there me discarded. The suc- 
cessive patches are merged smoothly within the band of 
agreement. The motivation for this approach is found in 
section 3 following a general discussion of the problem of 
planetary gravity and the information flow in our analysis, 
which are presented in section 2. 

2. Venus Gravity and Topography 

Figure 1 shows the natural flow of information from 
internal geophysical processes to our Doppler observable. 
Those internal processes, acting within the bulk of the planet 
and driven by the escape of heat, give rise to its mass 
distribution, which is partially manifested as topography. 
The mass distribution, including the topography, produces 
the external gravity field. Given the initial position and 
velocity of a spacecraft, its subsequent motion is determined 
by the Venus gravity field along with perturbations from the' 
Sun and other planets. The phase-coherent radio tracking 
system used by the NASA Deep Space Network yields a 
Doppler shift observable which reflects the spacecraft mo- 
tion. 

Ideally, we would like to invert all of these processes. 
Starting with the spacecraft observable, we would like to 
know the internal processes within the planet Venus. By 
well-established techniques, one can use the Doppler ob- 
servable to determine the spacecraft motion. Similarly, to 
within resolution limits imposed by the data, the observing 
geometry, and the noise, that same observable can be used 
to determine directly the external gravitational field. On the 
other hand, the internal mass distribution cannot, in princi- 
ple, be determined from the external gravitational field even 
if the latter were perfectly known. The problem of determin- 
ing the internal processes responsible for the mass distribu- 
tion is not even well defined. It comprises a major portion of 
the study of planetary interiors. 

To shed some light on those internal processes and their 
resulting mass distributions, we can construct models. For 
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Fig. 1. The flow of information (causality) from internal pro- 
cesses, which are driven by the temperature difference between the 
surface and the interior, to the spacecraft Doppler observable. Open 
arrows show mathematically possible inversions. 


those models applicable in the near-surface region, the 
observable consequences will include both the topography 
and the external gravitational field. For some types of 
internal processes, it is possible to obtain a relationship 
between the observed topography and the corresponding 
gravity anomalies. Such a signature can then be sought in the 
data. However, for Venus the resolution and fidelity of the 
gravity model are poor compared to that of the topography. 
Thus a direct comparison between the two would suffer a 
wavelength-dependent bias which would systematically dis- 
tort any conclusions about the internal structure. 

To overcome this problem of bias, we need a topographic 
map having the same distortion and spectral characteristics 
as the gravity model. Figure 2 is a simplified diagram of the 
information flow within our data processing system. In order 
to produce the required smoothed topography map, we start 
with the topographic data and calculate the corresponding 
spacecraft acceleration along the actual spacecraft trajec- 
tory, assuming the topography to be uncompensated and of 
unit density. These accelerations are processed by our linear 
inverter in the same way as the spacecraft Doppler rate 
residuals, yielding, respectively, the smoothed topography and 
gravity models. As shown by Goldberg and Reasenberg 
[1985], the resulting gravity model is of high fidelity; its 
location-dependent resolution and limited distortions are both 
well mimicked by the corresponding smoothed topography. 

3. Pioneer Venus and the Analysis Problem 

The orbital characteristics of PVO are the dominant factor 
in determining the resolution of our gravity maps and in 
selecting the analysis techniques. During the first 600 days of 
the mission, the altitude of the PVO at periapsis was 
maintained between approximately 140 and 190 km by 
means of propulsive maneuvers. The characteristics of the 
gravity maps that can be derived from the PVO tracking data 
are delimited by this low periapsis altitude, the high orbital 
eccentricity e = 0.84, an orbital inclination of 105° to the 
Venus equator, and the evolving Venus-Earth-Sun geome- 
try, which includes superior conjunctions and occultations 
by Venus of the spacecraft near periapsis. In particular, the 
PVO Doppler data support high-resolution maps in the 
vicinity of the spacecraft periapsis, whose latitude remained 
between 17° and I4°N during the first 600 days. However, at 
45° along track from periapsis. where the spacecraft local 


altitude was about 1 100 km, the possible resolution of the 
maps is correspondingly lower. Data used in our analysis 
were obtained when the Earth-Venus-spacecraft angle at 
periapsis was less than 90° and when the Earth-Sun-Venus 
angle was less than 150°. 

Two sets of circumstances together set the framework for 
selecting the analysis technique. The first pertains to the way 
the spacecraft moved, and the second pertains to where it 
moved with respect to the surface of Venus. These are 
addressed in order. The deterministic trajectories of celestial 
mechanics are an excellent representation of the motions of 
most planets and natural satellites on a time scale of re- 
corded history. These classical models justifiably neglect 
such stochastic perturbations as the variation in solar radi- 
ation pressure due to albedo variations. However, space- 
craft are different; they have about an eight-order larger 
area-to-mass ratio and may contain outgassing components 
(e.g., in the attitude control system) that produce accelera- 
tions significant to the analysis of tracking data. Such 
stochastic contributions to the driving terms of a model of a 
system are known generally as “process noise” and are 
dealt with by means of the Kalman-Bucy filter. (For a 
discussion of such filters, see, for example, Jazwinski [1970] 
and Gelb [1974].) 

The spacecraft acceleration along the line of sight is mea- 
sured by the Doppler rate data obtained by (numerical) differ- 
entiation of the Doppler tracking data. To be of use in a 
geophysical analysis, this measure must be “downward con- 
tinued” to the surface and transformed into a uniform repre- 
sentation. The surface resolution of such an inversion is limited 
to being not much better than A = 2trh, where A is the spatial 
wavelength on the surface and h is the local spacecraft altitude. 
For a surface mass density (represented as a surface harmonic 
series) the external vertical acceleration is 

av K r n+i " 

a= - =2 ( ' I+,) “TTJ 2 P «m (COS0) 

dr n~2 r m=0 

’ [ C nm cos ( m<p ) + S nm sin (m«p)] (1) 
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Fig. 2. The flow of information within our data processing 
system. The processors (computer programs) are in the boxes. The 
■•raw" data enter at the top. and the useful representations of the 
information are shown at the bottom as the products of the analysis. 
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where r, 6 , <p is the field point and r 0 is the reference radius 
taken to be the planet’s mean radius. The spatial wavelength 
corresponding to harmonics of degree n is A = 2irr 0 /n. For 
large n and for an altitude h, low compared to r 0 , 

(r/r 0 )- (n+2) = ^ -(n+2)A/r 0 ^ e —2irhl\ (2 ) 

The same result is easily obtained in the “flat planet” 
approximation. (For example, see Goldberg and Reasenberg 
[1985].) 

For the PVO analysis the available resolution varies 
substantially over the planet’s surface, principally as a 
function of latitude. The best resolution is at the latitude of 
periapsis: /i min > 140 km; 2ir(140 km) = 8.2° on the surface. 
To model the gravity at this resolution over the entire 
surface would require about 2.6 x 10 3 parameters, indepen- 
dent of the (nonredundant) form chosen for the gravity 
representation. Further, because of the variability of h, such 
a large number of parameters could not be independently 
determined. An attempt to do so would result in a degenerate 
(i.e. , severely ill-conditioned) estimator. 

Taking the apparently conservative approach of estimat- 
ing a model of uniformly low resolution is not an entirely 
satisfactory solution for two reasons. First, such a model 
would fail to represent all of the available information. 
Second, the unmodeled, high-spatial-frequency signatures in 
the gravity would be aliased into the modeled, lower- 
frequency signatures. Thus we are forced to the conclusion 
that the ideal representation would permit the resolution of 
the model to be variable over the planet surface such that it 
could reflect the intrinsic resolution of the data. 

The motions of Earth and Venus cause the observing 
geometry to vary with time. (Precession of the spacecraft 
orbit is relatively slow.) An analysis technique that requires 
a special observing geometry will, at best, be applicable to a 
small fraction of the available data and therefore should not 
be considered. However, even the most robust estimators 
show a decreased information rate for certain unfavorable 
geometry including an Earth-Sun- Venus angle near 180°. 

The classical approach of the celestial mechanician is to 
expand the central body potential as a spherical harmonic 
series and to analytically determine the effects of the indi- 
vidual terms of the series on the spacecraft orbital elements. 
The spacecraft orbit would then be determined from each of 
several short spans of data and the evolution of the elements 
used to estimate the harmonic series coefficients. For the ab 
initio analysis of the PVO tracking data, this approach has 
two fatal flaws. First, the spherical harmonic representation 
yields a uniform surface resolution which, as previously 
discussed, is inappropriate to the PVO problem. Second, the 
use of orbital elements discards too much of the information 
content of the data. 

A Venus gravity model was determined by Ananda et al. 
[1980] from the analysis of the evolving spacecraft orbital 
elements. A principal motivation for this work was the 
availability of the spacecraft elements as determined each 
day by the PVO Navigation Team. The cost of the final 
analysis was thus relatively low, making this an attractive 
type of study. However, the model is of low resolution. 

Early planetary and lunar orbiter gravity analyses in- 
cluded the direct, least squares estimation of harmonic 
coefficients from the Doppler tracking data. (Attempts to use 
ranging data were generally not successful.) It was widely 
believed that the best results would be obtained from the 


analysis of the longest possible continuous spans of data. 
This belief appears to have been based on an examination of 
the variational equations which showed increased sensitivity 
with time. 

In determining the gravity field of Mars from the Mariner 
9 tracking data, Reasenberg et al. [1975, p. 89] showed that 
the direct, least squares analysis yielded better (i.e., more 
consistent) results when the data spans were broken, even 
though this required increasing the parameter set to include 
additional sets of spacecraft elements. The work was moti- 
vated by a concern for the effects “of important random, or 
quasi-random, accelerations of the spacecraft, due, for ex- 
ample, to imbalances in the gas jets used to control the 
orientation of the spacecraft . . The multiple short-arc 
analysis was viewed as an approximation to a Kalman filter 
which, at the time, had neither been implemented in our 
software, nor used by any group to- estimate planetary (or 
lunar) gravity. 

For a Venus gravity model based on the PVO tracking 
data and spherical harmonics, the previously discussed 
conflict over resolution would be particularly severe. A 
possible solution to this conflict would be to introduce a 
priori constraints on the estimated parameters such that the 
effective resolution would vary over the surface in corre- 
spondence to the intrinsic resolution of the data. Although 
this is possible in principle, we know of no case in which this 
approach has been applied to the determination of a gravity 
model. Even if it were to be applied, in the case of PVO, it 
would result in the need for computationally burdensome 
data processing. In the absence of such a constraint, the 
model would likely contain artifacts not usefully related to 
the actual gravity of the planet. For example, these might 
appear as features of small lateral extent in a region where 
the data are unable to support such resolution, i.e., away 
from the periapsis latitude. The alternative of using a diag- 
onal a priori constraint yields an undesirable bias to all 
harmonic coefficients. 

The earliest use of short-arc analysis was by Muller and 
Sjogren [1968]. They mapped the Doppler residual rate 
directly to the surface of the Moon along the line of sight. 
This direct residual mapping (DRM) technique permitted 
them to make the first detection of the lunar mascons. DRM 
has been used extensively by W. L. Sjogren and his cowork- 
ers [e.g., Sjogren et al., 1974, and references therein]. The 
technique has proved of great value for a first, qualitative 
examination of gravity, since it yields maps of good resolu- 
tion with a minimum of data processing. However, because 
there is no inversion involved, the resulting gravity maps are 
ill-suited for direct quantitative interpretation as noted (often 
and with clarity) by Sjogren, who cites numerical studies by 
Gottlieb [1970]. 

Of course, in the same way one could use the Doppler 
data, it is possible to use a DRM to fit a specific geophysical 
model, as was done by Phillips et al. [1978]. Thus one might 
be persuaded that all gravity representations are equally 
good since to correctly test a specific geophysical model 
requires calculating the corresponding external gravity map 
that would have been determined by the method used for the 
available gravity model. However, an inversion which main- 
tains the full resolution of which the data are capable 
provides better discrimination among geophysical models; a 
clean inversion better supports the synthesis of useful hy- 
potheses and their preliminary testing. 
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Over the years, several ad hoc schemes for inverting 
Doppler gravity data have been proposed. (For example, see 
Bowin [1983]. In his Figure 9 the two sets of discontinuities, 
which follow spacecraft tracks and are each due to a 
transition from one set of contiguous tracks to another, 
attest to the failure of the technique.) Such schemes may 
yield maps that give the appearance of being sensible but 
which can make no contribution to the knowledge of the 
internal structure of the body. (To the extent that they are 
used for further analysis, such maps may interfere with 
progress.) This is particularly regretable since the literature 
contains numerous good papers on the geophysical inverse 
problem including Backus and Gilbert [1970], Burkhard and 
Jackson [1976], Moritz [1976, 1978], and Jackson [1979]. 

If the ad hoc schemes are incorrect, why do the maps look 
plausible? An inversion is generally a linear operation on the 
data to yield a gravity map. The Doppler tracking data are 
“closely related” to the gravity as evidenced by the success 
of the DRM. A plausibly useful linear operation on a DRM 
would result in location-dependent shifts in the amplitude 
and phase of the components of the map but would not be 
likely to obscure or reveal major features. The same is true 
of a linear operation on the Doppler residuals. It is hard to 
imagine an inversion technique being proposed and found to 
obscure the major features of the gravity field. However, 
there is a considerable difference between not obscuring 
features and yielding a map suited to quantitative interpre- 
tation. The former is hardly an achievement; the latter was 
the objective of our analysis. 

4. The Kalman Filter in PEP 

The Kalman filter (KFj in PEP has a few unusual features. 
These and some other pertinent aspects are described here 
functionally; the details of the implementation are driven by 
the peculiarities of PEP. In section 5, we present the results 
of several numerical tests of the PEP KF. 

In the classical KF, the state estimate is updated with each 
datum. The state estimate and covariance are propagated 
forward, and the effect of the process noise is included in 
preparation for the next datum. In this way, the optimal 
estimates are available shortly after the receipt of each 
datum, a considerable advantage for real-time analysis and 
control. (For a recent astronomical application of a real-time 
filter, see Reasenberg [1990].) However, for the analysis of 
scientific data, immediacy is not important. Therefore, in the 
PEP KF, the data are collected into batches and each batch 
is used to form a set of normal equations as would be done 
for a weighted-least-squares analysis. If these normal equa- 
tions were added together directly, we could obtain the 
standard weighted-least-squares solution. 

The batch normal equations are processed in place of 
individual data in the PEP KF. The state is defined such that 
the nominal state propagation matrix is the identity. Thus, 
instead of including the position and velocity of a spacecraft 
in the state, the corresponding initial conditions are in- 
cluded. Since the process noise covariance is defined in 
terms of the spacecraft position and velocity as discussed 
below, the scheme requires that a set of variational equa- 
tions be integrated and that these be used to map the process 
noise to the initial epoch. (In the more common form of the 
KF, those same variational equations would be required to 
calculate the state transition matrix.) 


In the PEP KF analysis of the PVO tracking data, we have 
an advantage not usual in the KF analysis. The PVO 
Navigation Team used single-day batches of data to deter- 
mine the spacecraft state for most of the first 600 days of the 
orbital phase of the mission. Although these state estimates 
were found to differ slightly from those derived by PEP, as 
initial estimates, they proved indispensable to the efficient 
conduct of our analysis. The PEP KF contains a provision to 
accept initial state estimates at several epochs during the 
time span under consideration. In the first step of the 
analysis, these are used to find a reference trajectory. 
Starting with the earliest epoch, the trajectory is numerically 
integrated until the next epoch. Here an integration transi- 
tion is performed: The difference between the integrated and 
externally supplied states is found and mapped to the initial 
epoch using the variational equations. The current state 
difference is used to modify the difference tables of the 
numerical integration so that the trajectory passes through 
the new externally supplied point and continues. The state 
offset is mapped to the initial epoch and saved to be used 
with the linearized estimator. This process is repeated for 
each of the state vectors until the requested end of the 
integration is reached. 

The atmospheric drag at periapsis makes the PVO trajec- 
tory analysis significantly less linear than the corresponding 
drag-free problem. In our early analyses, before the imple- 
mentation of the multistate starting procedure, the trajectory 
fitting typically required six iterations per batch of data. 
With the multistate starting procedure, the initial (segment- 
ed) trajectory proved to be so close to “correct” that no 
iteration was necessary as will be discussed in the next 
section. Doppler residuals, suitable for geophysical inver- 
sion, were obtained by linear prediction from the prefit 
residuals, the state adjustment vectors of the first KF 
solution, and the sensitivity matrix (i.e., dzlda, where z is 
the observable and a is the vector of parameters). 

The process noise model consisted of two parts. The first 
applied to the atmospheric density parameter and was made 
large enough that the density estimates for adjacent days 
were essentially independent. The second part of the model 
was an isotropic acceleration variance density. Each Carte- 
sian component had an amplitude of 10“ 19 AU 2 /d 3 . 

5. Kalman Filter Performance 

We have conducted two types of tests of the PEP KF as 
used for the analysis of the PVO Doppler data. In the first set 
of tests, we used Doppler residuals kindly provided to us by 
W. L. Sjogren of the Jet Propulsion Laboratory (JPL). We 
compared these to the corresponding PEP KF residuals. In 
the second set of tests, PEP was used to numerically 
integrate a PVO trajectory with one of several models and to 
calculate the expected Doppler data. The PEP KF was used 
to analyze these simulated data and the resulting postfit 
residuals, or derived Doppler rate residuals, were investi- 
gated. Below we discuss the two sets of tests. 

JPL Comparison 

Sjogren provided us with residuals that he had obtained 
from his own (PVO) spacecraft orbit determinations. In this 
work, he used 2-hour spans of tracking data starting 1 hour 
before periapsis and ending I hour after periapsis. The 
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trajectory was integrated using an atmospheric model deter- 
mined for each orbit by the Navigation Team. 

A comparison of the Doppler residuals showed a system- 
atic difference that is no larger than the data noise. This 
difference has no apparent high-frequency component and 
thus no effect on the gravity model. (Note that the data noise 
cancels in this comparison.) We concluded that there would 
be little difference in our maps were we to replace our KF 
residuals with Sjogren’s short-arc batch-analysis residuals. 

Simulation Study 

We have done a series of numerical experiments that 
address the accuracy of the PEP KF analysis of the PVO 
data. In each experiment, PEP was used first to integrate a 
spacecraft trajectory with a given set of parameters and 
second to calculate the corresponding Doppler observables. 
The PEP KF was used to estimate a ‘‘best fit” filter 
trajectory assuming no anomalous gravity and a nominal 
atmospheric scale height. The resulting Doppler residuals 
were numerically differentiated to determine the Doppler 
rate residuals, which were expected to be proportional to the 
“line-of-sight” (LOS) component of the unmodeled part of 
the spacecraft acceleration. 

In each study discussed below, we simulated the same set 
of seven contiguous orbits (beginning and ending at 
apoapse), with initial conditions (IC) and epoch taken from 
an actual PVO orbit. The times and IC for each KF epoch 
were obtained from the simulations, with IC rounded to from 
7 to 10 figures, to mimic the typical accuracy of the (JPL) IC 
used in actual PVO/KF orbit determinations (OD). In actual 
OD, perturbations due to solar gravitation and solar radia- 
tion are modeled. In both the simulations and the OD for the 
numerical studies, the former perturbation is included, the 
latter is not. Planetary perturbations of the spacecraft, which 
are not included in the analyses of the PVO tracking data, 
were not included in the simulations. 

In the first numerical experiment, the initial integration 
assumed neither atmospheric drag nor anomalous gravity. 
Thus the KF deterministic model was correct. The resulting 
LOS residuals were systematic (presumably due to imper- 
fect adjustment of the rounded IC) and had a RMS of under 
0.0003 mGal and a peak of 0.0012 mGal in the periapsis 
region. This is about 4 orders of magnitude below either the 
typical error due to the inversion or the worst errors encoun- 
tered in the other KF tests. 

In the second numerical experiment, the initial integration 
model included three gravity-harmonic terms (7 3 = 5 x 
10 -5 , 7| 0 = — 3 x 10 -5 , 0 * 20,5 = 1 x 10 —5 ) and zero 
atmospheric density. The LOS residual is shown in Figure 
3 a along with the difference between this and the expected 
LOS acceleration calculated directly from the gravity- 
harmonic terms. This difference, which is on a 10X finer 
scale, does not appear to include the signature of the gravity 
model; the resulting error is presumed to be under 1%. 

In the third numerical experiment, the gravity-harmonic 
terms were zero but atmospheric drag was included in the 
initial integration. The atmospheric density p at altitude z 
was modeled as 

p(z) = Po^" U “' o)/ " (3) 

where p 0 is the density at the reference altitude z 0 and H is 
the scale height. For the integration we used H = 7 km, z 0 



Fig. 3. Results of tests of the performance of the Kalman filter 
(KF). These test address the question: To what extent does the filter 
yield the line-of-sight (LOS) acceleration of the spacecraft? See text 
for a description of the tests. Solid curve shows LOS residual; 
dashed curve shows difference between LOS residua) and expected 
acceleration due to gravity harmonics, (a) Experiment two, recov- 
ery of gravity harmonics, (b) Experiment three, effect of atmo- 
spheric drag, (c) Experiment four, combination of gravity harmonics 
and atmospheric drag. 


= 150 km, and p 0 = 5 x 10 -13 g/cm 3 . In the KF analysis, we 
used a scale height of H' = 10 km and the previous value for 
z 0 . We selected Po such that for the first periapsis passage, 
p(z p ), the density at periapsis, yielded about the correct 
value of p p H il2 , since this value would tend to keep the 
problem from becoming excessively nonlinear. (A similar 
approach was used with the real data; the value of piz p )H' n 
was obtained from the orbital elements determined by the 
JPL Navigation Team.) 

Figure 3 b shows the LOS residual from the KF fit. The 
signature is that expected from the use of a “wrong” scale 
height in the KF. Table 1 (experiment 3) shows the determi- 
nation of the atmospheric density by the KF. It is easily 
shown that the change of spacecraft orbital period due to 
drag is proportional to y = p(z p )H V2 . Table 1 shows that 
this quantity is estimated with a small fractional error 
although the error is large compared to the standard devia- 
tion <7. Based on the analysis of the PVO data, we know that 
the density shows much larger fluctuations. 

The mismatch in scale height used in this test is larger than 
we expect most of the corresponding mismatches to be in the 
actual data analysis. Further, in the geophysical inversion, 
information from several spacecraft passes is combined 
(averaged) to produce the gravity estimate. Thus unless the 
scale height is systematically incorrect for several days, 
there should be some reduction of the effect of the atmo- 
sphere mismodeling on the gravity model. We therefore 
expect the atmospheric modeling errors to result in gravity 
model errors of well under 5 mGal. 


ERROR IN FILTER RESIDUAL (mgal) 
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TABLE 1. Determination of Atmospheric Density by the Kalman Filter 


Orbit 

Z P z 0> 

km 

Zp Zq » 

km 

P°. 

10" 13 

g/cm 3 

Piz p ), 

10" 13 

g/cm 3 

P(z p ), 

10~ 13 

g/cm 3 

>’ = P/ z p) 

H 1/2 

y-y 

a 

(y - y)/cr 

1 

0.15 

0.15 

4.12 

Experiment 3 
4.89 4.06 

12.9 

-0.114 

0.0064 

-18 

2 

0.79 

0.79 

4.01 

4.46 

3.71 

11.8 

-0.090 

0.0063 

-14 

3 

1.57 

1.57 

3.87 

4.00 

3.31 

10.6 

-0.117 

0.0062 

-19 

4 

2.48 

2.47 

3.71 

3.51 

2.90 

9.3 

-0.117 

0.0060 

-20 

5 

3.53 

3.53 

3.53 

3.02 

2.48 

8.0 

-0.148 

0.0059 

-25 

6 

4.72 

4.72 

3.37 

2.55 

2.10 

6.7 

-0.100 

0.0059 

-17 

7 

6.06 

6.06 

3.17 

2.10 

1.73 

5.6 

-0.106 

0.0060 

-18 

1 

0.25 

0.27 

4.31 

Experiment 4 
4.82 4.19 

12.8 

0.502 

0.0065 

78 

2 

1.01 

1.05 

4.20 

4.33 

3.79 

11.4 

0.527 

0.0063 

84 

3 

1.91 

1.96 

4.06 

3.81 

3.34 

10.1 

0.478 

0.0062 

77 

4 

2.94 

2.98 

3.87 

3.28 

2.88 

8.7 

0.407 

0.0060 

68 

5 

4.12 

4.14 

3.71 

2.78 

2.45 

7.3 

0.398 

0.0060 

67 

6 

5.43 

5.45 

3.55 

2.30 

2.06 

6.1 

0.412 

0.0059 

70 

7 

6.89 

6.93 

3.40 

1.87 

1.70 

4.9 

0.432 

0.0060 

72 


The model of the atmospheric density p and its parameters p 0 = 5 x 10 -13 g/cm 3 ,z 0 = 150 km, and H = 7 km, are defined in the text. 
The variables with a circumflex (e.g., p) refer to quantities estimated (directly or indirectly) in the KF orbit determination; the corresponding 
unmarked variables refer to their “true” values, i.e, those used in or obtained from the original simulation. z„ is the spacecraft altitude at 
periapse, y is the quantity proportional to the change in orbital period, H = 7 km and H' = 10 km are the scale heights used in the simulation 
and KF, respectively, and cr is the formal error for y . 


In the fourth numerical experiment, the initial integration 
combined the gravity model of experiment 2 with the atmo- 
spheric model of experiment 3. The KF again used the 
“wrong” scale height, H' = 10 km. Figure 3c shows the 
LOS residual and, again on a 10X finer scale, the difference 
between this residual and the expected contribution due to 
the harmonics. Table 1 (experiment 4) shows the density 
determination by the KF. Although the errors in the density 
estimates (see y-y in Table 1) are 3-6 times those of the 
previous case, they are still small compared to the day-to- 
day fluctuations seen in the PVO-derived density estimates. 
Thus it is better to estimate p for each periapsis pass than to 
assume po constant, even for a few days. 

A comparison of the error in the KF accelerations of 
Figure 3c with those of Figures 3a and 3 b shows that the 
error produced by a combination of unmodeled anomalous 
gravity and incorrectly modeled atmospheric density is, to a 


good approximation, simply the sum of the errors due to 
each model error individually. This, as well as the results 
presented in Table 1 , indicates that the low-frequency grav- 
itational harmonics are not significantly “absorbed” or 
“aliased” as atmospheric process noise. 

We extended the above numerical experiments to examine 
the errors introduced by not iterating the estimator. In the 
iteration associated with experiment 2 (gravity harmonic 
terms but no atmospheric drag), the estimator required two 
iterations to converge plus one to confirm convergence. We 
found in experiment 3 (drag but no harmonics) that, as 
suggested by previous work, the atmospheric drag term made 
the problem highly nonlinear. (This reflects the small scale 
height of the atmosphere.) The same nonlinearity was present 
when the gravity harmonic terms were added (experiment 4); 
the gravity terms caused larger adjustments to the elements 
and thus tended to hide the nonlinearity. The iteration showed 



EAST LONGITUDE 

Fig. 4. Spacecraft tracks projected on the surface of Venus and shown 5° past the upper and lower edges of the 
maps in the color plate. For the gravity inversion, data were used from somewhat beyond the region shown here: in 
particular, we excluded all data taken with a spacecraft altitude of more than 4000 km. 






wmmm 


srutira 


Reasenberg and Goldberg: High-Resolution Gravity Model of Venus 


Plate 1. The gravity map (above) and corresponding smoothed topography map (below) of Venus. These maps 
cover.2/3 of the surface of Venus and are cut off where the resolution has become poor. The false color of the maps 
correspond to altitude. For the gravity map, it is the altitude to which material of density half the mean density of Venus 
would need to be piled to yield the observed external potential. The color bar between the maps shows increasing 
altitude from left to right in steps of 150 m and 333 m for the gravity and topography maps, respectively. Near periapse, 
the discretization limits the resolution of the gravity map to A = 4°. However, the combination of the exponential loss 
of signal (equation (2)) and the behavior of the estimator cause the response to roll off as the wavelength gets small. We 
estimate the effective resolution to be A = 6° at periapsis and to be larger by a factor of 10 at the upper and lower edges 
of the map. By construction, the topography map should have the same resolution characteristics (and the same 
distortions) as the gravity map. 


that the Doppler rate residuals from the converged solution 
differed systematically from the ones predicted linearly from 
the first iteration by about 0.2 mGal RMS; no damage is done 
to our gravity maps by not iterating. 


6. Data Analysis and Results 


Plate 1 shows false-color images of our gravity map (top) 
and the corresponding smoothed topography map (bottom). 
These have been scaled to provide easy comparison. The 
maps are based on the combined results of nine separate 
inversions, and incorporate 1.2 x 10 5 Doppler rate data from 
251 spacecraft revolutions, which occurred between April 
1979 and August 1980. Displayed in Figure 4 are projections 
on the Venus surface of the spacecraft trajectory at the 
transpond times of the Doppler data we used. Note that the 
gravity map shows no evidence of either the data gaps or the 
redundancy near45°E longitude. The orbits were determined 
in 38 separate Kalman filter batches, generally comprising 
six to eight orbital revolutions each. The breaks between 
batches were usually dictated by the occurrence of propul- 
sive spacecraft maneuvers. 

We model the external gravity as the sum of a point mass 
centered on the planet and a surface mass density. In each 
inversion region, we have discretized the planetary surface 
with lines of constant latitude or longitude. This gives rise to 
trapezoidal cells (neglecting curvature), which we take to 


TOPOGRAPHY 


Fig. 5. Scatter of gravity (vertical axis) versus topography, from 
the same data displayed in Plate 1 . taken in the band of highest 
resolution, between 5°S and 35°N over the full longitude range at 1° 
lattice points. The uniis for the respectiv e quantities are the same as 
in Plate I . 





Scatter plots of gravity and topography, as in Figure 5, but front 40° x 40° regions; (a) 0°-40°E; (b) 40°-80°E; (c) 80°-120°E; (d) 
120°-160°E; (e) 16CT_200°E; (/) 200°-240°E; (g) 240°-280°E; (h) 280°-320°E; («) 320°-360°E. 
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have uniform surface mass densities. At periapsis, the cells 
are 2° x 2°, and for the highly eccentric orbit of PVO, we 
advantageously make the linear dimension of each cell 
roughly proportional to the local spacecraft altitude. The cell 
is divided into four (i.e., 2x2) “subcells, ’’ and the field of 
each subcell is modeled as the sum of a centered point mass 
and quadrupole; the dipole term is neglected. 

As discussed by Goldberg and Reasenberg [1985], the 
data inversion is performed by a linear Bayesian least 
squares estimator in which we use the weighted Jekeli 
covariance to form an a priori covariance matrix. We set a 
500-parameter limit to the inversions, which permits a region 
covering about 60° of longitude. We established inversion 
regions that have a considerable overlap. Within each over- 
lapping portion is a strip in which there is not a significant 
difference between the two adjacent inversion regions. The 
data displayed at each point in the maps of Plate 1 either are 
taken directly from a single inversion region or are the 
weighted average of the data from two regions within such a 
strip of near agreement. In the latter case, the weighting 
varies linearly across the strip. 

A cursory examination of Plate 1 reveals what has long 
been known, that gravity is well correlated with topography 
at long spatial wavelengths. More careful examination re- 
veals that this correlation continues down to the smallest 
wavelengths resolved: virtually every visible feature in the 
smoothed topographic map has a gravitational counterpart. 
The reverse is not always true. However, many apparently 
unmatched small gravitational features do have counterparts 
in the unsmoothed topography [e.g., Pettengill et al., 1980]. 
For example, the principal peaks and central depression 
seen in the gravity over Asteria Regio (around 20 8 N, 265°E) 
all correspond to real topographic features that are not seen 
in the smoothed topography map. A similarly strong exam- 
ple can be seen in the detailed structure of Atla Regio 
(5°-25°N, 180°-200°E). 

It appears that most gravitational features that are at least 
two contour levels high have a real topographic counterpart. 
The preferential absence of many of these counterparts from 
the filtered topography for short spatial wavelengths is an 
indication that the spectral admittance on Venus tends to be 
higher at shorter wavelengths than at longer ones [cf. Gold- 
berg and Reasenberg, 1985]. That many single-contour-level 
gravitational features do not have counterparts in the topo- 
graphic data is consistent with our estimate of the noise level 
of the gravity map which is about one contour. 

It has been apparent from the low-frequency gravity 
studies that the level of isostatic compensation on Venus is 
far from uniform. In fact, the term compensation may be 
misdescriptive since at least some of the features may be the 
result of dynamic processes [Kiefer et r '. . 1986]. Some large 
topographic features of unusually high elevation, such as 
Beta Regio (25°N, 280°E) and Ulfram Regio (5°N, 200°E), 
correspond to similarly extreme gravitational features, while 
other high topographic regions, such as those in western 
Aphrodite Terra (5°S, 85°E and 10°S, 135°E), correspond to 
relatively modest gravity peaks. In Plate 1, such contrasts 
can also be seen for much smaller features, such as the 
components of Phoebe Regio (0°-20°S, 275°-300°E), Eisila 
Regio (10°-30°N, 10°W-30°E), and the region west of Beta 
Regio (0°— 45°N. 210°-260°E). 

There appear to be two distinct trends of isostatic com- 
pensation. as evidenced by the buurcated relation between 


gravity and topography as seen in Figure 5. In the units of 
our maps, one trend in this scatter plot has a gravity-to- 
topography ratio of about 1/3; the other, about 1/10. This 
bimodality occurs locally as well as planet wide. Figure 6 
shows the slopes of 1/3 and 1/10 appearing in smaller regions, 
with little or no evidence for other slopes. Of particular 
interest is the appearance of the 1/3 slope in Figures 6a and 
6b, where it occurs in regions between 0 and 2 km elevation. 
Such “rolling plains” regions have been thought to be more 
fully compensated than other portions of the planet, but this 
appears not always to be the case. 

7. Summary 

We have developed techniques that permit the efficient 
and accurate determination of the gravity field of a planet 
from the Doppler tracking of an orbiting spacecraft. These 
techniques work in the case of a highly eccentric orbiter and 
permit a gravity model whose resolution varies over the 
surface of the planet as required, for example, corresponding 
to the local spacecraft altitude. 

We have applied these techniques to the Doppler tracking 
data from the Pioneer Venus Orbiter to yield a Venus gravity 
model. The gravity map appears devoid of the computational 
artifacts that have traditionally plagued planetary gravity 
inversions. In particular, all of the features of the gravity 
map that are larger than the noise level correspond to 
topographic features. Since no topographic information was 
used in the inversion of the Doppler data to yield the gravity 
map, the strong correspondence is taken to be real and to 
indicate a lack of artifacts. (The conspiracy theory is reject- 
ed.) 

In order to permit an easy comparison of the gravity and 
topography, we have generated a topography map with the 
same resolution and distortion as the gravity map. This 
smoothed topography map was made by (1) converting the 
high-resolution topographic data from the PVO altimeter 
into pseudo-Doppler data that correspond to the real Dop- 
pler data and (2) performing the same inversion on these as 
on the real Doppler data. Scatter plots of gravity against 
smoothed topography show two distinct trends, suggesting 
two modes of topographic support or history. More detailed 
analysis of this phenomenon, for example by means of 
spectral admittance, seems warranted. 

Note that the digital data that were used to make the two 
maps in Plate 1 have been provided to the Magellan Project 
for their preencounter data set and to the NSSDC for 
distribution. 
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